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ABSTRACT 

While many geological and geophysical processes such as the melting of icecaps, the magnetic expression of 
bodies emplaced in the Earth's crust, or the surface displacement remaining after large earthquakes are spatially 
localized, many of these naturally admit spectral representations, or they may need to be extracted from data 
collected globally, e.g. by satellites that circumnavigate the Earth. Wavelets are often used to study such 
nonstationary processes. On the sphere, however, many of the known constructions are somewhat limited. And 
in particular, the notion of 'dilation' is hard to reconcile with the concept of a geological region with fixed 
boundaries being responsible for generating the signals to be analyzed. Here, we build on our previous work on 
localized spherical analysis using an approach that is firmly rooted in spherical harmonics. We construct, by 
quadratic optimization, a set of bandlimited functions that have the majority of their energy concentrated in an 
arbitrary subdomain of the unit sphere. The 'spherical Slepian basis' that results provides a convenient way for 
the analysis and representation of geophysical signals, as we show by example. We highlight the connections to 
sparsity by showing that many geophysical processes are sparse in the Slepian basis. 

Keywords: spectral analysis, spherical harmonics, statistical methods, geodesy, inverse theory, satellite geodesy, 
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1. THE SPHERICAL SLEPIAN BASIS 

We denote the colatitude of a geographical point f on the unit sphere surface r^ = {f: ||f|| = l}byO<6><7r and 
the longitude by < ^ < 27r. We use R to denote a region of 1^, of area A, within which we seek to concentrate 
a bandlimited function of position f = (O^cj)). We use orthonormalized real surface spherical harmonics,^' ^ thus 
expressing a square-integrable real function /(f) on the surface of the unit sphere as 
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/(f) = /^m>lm(r), flm = / fYimdVt, and / YimXi'm' dVt = 5w5mm'- (1) 
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The Slepian basis for the domain R is the collection of bandlimited functions 

L I 



(2) 



^(f ) = gimyimi^) for which ^= 9^{'^)dn I (f) = maximum. 

Maximizing equation ([2| leads to the spectral-domain Hermitian, positive-definite eigenvalue equation 

L I' . 
^ ^ ^Im.Vm'gVm' Mlm, with Dim.l'm' \ yimXv m' dVt, < I < L, (3) 

7/ n 7/ ^ R 



l'=0 m' = -l 

but we may equally well rewrite eq. ([s]) as a spatial-domain eigenvalue equation: 

[ D{r,r')9{^')dn' = Xg{r), with Z)(f,f') = E (^)^'(^ ' f')' ^ ^ (4) 

where Pi is the Legendre function of integer degree /, which arises in this setting as a consequence of the spherical 
harmonic addition theorem. Eq. Q is a homogeneous Fredholm integral equation of the second kind, with 



a finite-rank, symmetric, Hermitian kernel, and tiie finite set of bandlimited spatial "Slepian" eigensolutions 
^i(r),^2(r), . . . ,^(L+i)2(f) is orthonormal over the whole sphere Q and orthogonal over the region R: 



/ gagpdQ = Saf3, and / gagf3 dQ = XaSafd- (5) 
Jn JR 

When the concentration region is a circularly symmetric cap of radius centered on the North Pole the solutions 
to eq. Q break down by order m and are separable in 6 and ^, and the colatitudinal parts g{9)^ which depend 
only on |m|, are identical to those of a Sturm-Liouville equation which can be solved in the spectral domain by 
diagonalization of a simple tridiagonal matrix with an almost linear spectrum.^' ^ We define a space-bandwidth 
product or 'spherical Shannon number' by the sum of the eigenvalues, 

N= ^« = E E Dim,im= / D{r,r)dn = {L^lf^. (6) 



a=l 1=0 m= — l 



The complete set of bandlimited spatial Slepian eigenfunctions ^1,^2, • ,^(l+i)2, irrespective of the particular 
region of concentration that they were designed for, are a basis for bandlimited scalar processes anywhere on 
the surface of the unit sphere.^' ^ This follows directly from the fact that the spectral localization kernel ([s]) 
is real, symmetric, and positive definite: its eigenvectors giim^92im^ • • • ^9{L-\-iyim form an orthogonal set, thus 
the Slepian basis functions ^^(f), a = 1, . . . , (L + 1)^ given by eq. ^ simply transform the same-sized limited 
set of spherical harmonics yi^(f), 0<1<L, —l<m<l that are a basis for the same space of bandlimited 
spherical functions with no power above the bandwidth L. After sorting the eigenvalues in decreasing order, this 
transformation orders the resulting basis set such that the energy of the first N functions, ^i(f), . . . ,^Ar(f), with 
eigenvalues A ~ 1, is concentrated in the region i?, whereas the remaining eigenfunctions, ^Ar+i(f), . . . , ^(l+i)2 (f), 
are concentrated in the complimentary region R = Q — R. As in the one- and two-dimensional case,^^ therefore, 
the reduced set of basis functions ^1, ^2, • • • , ^at can be regarded as a sparse, global, basis suitable to approximate 
bandlimited processes that are primarily localized to the region R. The dimensionality reduction is dependent 
on the fractional area of the region of interest, i.e. the full dimension of the space {L-\- 1)^ can be "sparsified" to 
an effective dimension of N = {L -\- 1)'^A/{4:7t) when the signal of interest lies in a particular geographic region. 

An example of Slepian functions on a circular region on the surface of the sphere can be found in Figure [l] 

2. PROBLEMS IN GEOPHYSICS (AND BEYOND) 

With all of the foregoing established as fact and referring again to the literature cited so far for proof and further 
context, we return to considerations closer to home, namely the estimation of geophysical signals from noisy 
and incomplete observations collected at or above the surface of the spheres "Earth" or "planet". We restrict 
ourselves to real- valued scalar measurements, contaminated by uncorrelated additive noise, which we may not 
know but which we shall describe by idealized models. We focus exclusively on data acquired and solutions 
expressed on the unit sphere. We have considered generalizations to problems involving satellite data collected 
at an altitude and/or potential fields elsewhere.^' ^' ^' Two different statistical problems arise in this context, 
namely, (i) how to find the "best" estimate of the signal given the data,^ and (ii) how to construct from the 
data the "best" estimate of the power spectral density of the signal in question. In this contribution we limit 
ourselves to problem (i) as it is here that the connections to sparsity are most readily apparent. 

Thus, let there be some data distributed on the unit sphere, consisting of 'signal', n and 'noise', s, and let 
there be some region of interest R C ft^ in other words, let 

M^) ^ / ^(r) + ri{r) ifreR 
^ ' \ unknown/undesired if r ^ ft — R. ^ ' 

We assume that the signal of interest can be expressed by way of spherical harmonic expansion as in eq. ([T]) , and 
that it is, itself, a realization of a zero-mean, Gaussian, isotropic, random process, namely 

00 I ^ 

5(f) = ^ ^ SimYim{r), Sim = / sYimdQ, (sim) = and {simSUm') = Si6w6mm'' (8) 

^=0 m=-l 



For convenience we furthermore assume that the noise is a zero-mean stochastic process with an isotropic power 
spectrum, i.e. (n(f)) = and (n^^n//^/) = Ni Su^Smm'^ and that it is statistically uncorrelated with the signal. 
We refer to power as white when Si = S or Ni = or, equivalently, (n(f)n(f')) = NS{t,t^). Our objective is 
to determine the best estimate sim of the spherical harmonic expansion coefficients sim of the signal. While in 
the real world there can be no limit on bandwidth, practical restrictions force any and all of our estimates to be 
bandlimited to some maximum spherical harmonic degree L, thus of necessity sim = and Si = when I > L: 

L I 

^(f) = E E ^InrYlmi^). (9) 

^=0 m=-l 

This limitation, combined with the statements eq. ^ on data coverage and the region of interest, naturally puts 
us back in the realm of 'spatiospectral concentration'. As we shall see, solving the problem at hand will gain 
from involving 'localized' Slepian functions rather than, or in addition to, the 'global' spherical harmonics basis. 

This leaves us to clarify what we understand by "best" in this context. While we adopt the traditional 
statistical metrics of bias, variance, and mean squared error to appraise the quality of our solutions, the 
resulting connections to sparsity will be real and immediate, owing to the Slepian functions being naturally 
instrumental in constructing efficient, consistent and/or unbiased estimates of sim and/or Si. Thus, we define 

v = {s^)-{sf, b={s)-s, e = s-s, and (e^) = v + 6^ (10) 

where the lack of subscript indicates that we can study variance, bias and mean squared error of the estimate of 
the coefficients sim but also of their spatial expansion 5(f), or indeed of their power spectrum Si. 

Signal estimation from noisy and incomplete spherical data 
Spherical harmonic solution 

Paraphrasing results elaborated elsewhere,^ we write the bandlimited solution to the damped inverse problem 

/ {s — d)'^ dft -\- T] / (il] = minimum, (11) 
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where 77 > is a damping parameter, by straightforward algebraic manipulation, as 

L V „ 
Slm = ^ ^ {Dlm,l'm' ^VDlm,l'm')~^ dYi>m'dVt, (12) 

where Dira^Um'^ the kernel that localizes to the region R = Q — compliments Dim,Um' given by eq. ([s]) which 
localizes to R. Given the eigenvalue spectrum of the latter, its inversion is inherently unstable, thus eq. (11) is 
an ill-conditioned inverse problem unless 77 > 0, as has been well known, e.g. in geodesy. As can be easily 
shown, without damping the estimate is unbiased but effectively incomputable; the introduction of the damping 



term stabilizes the solution at the cost of added bias. And of course when R = eq. (12) is simply the spherical 
harmonic transform, as in that case, eq. ([3| reduces to eq. ([T]), in other words, then Dim,Um' = ^w^ram'' 

Slepian basis solution 

We could seek a trial solution in the Slepian basis designed for this region of interest R by writing 

(L+l)2 

'^(^) = X] ^c,^c.(r). (13) 

This would be completely equivalent to the expression in eq. ([9| by virtue of the completeness of the Slepian basis 
for bandlimited functions everywhere on the sphere and the unitarity of the transform ([2| from the spherical- 
harmonic to the Slepian basis. The solution to the undamped {rj = 0) version of eq. (ITT]) would then be 



dgadQ, (14) 
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which, being completely equivalent to eq. (12) for r] = 0, would be computable, and biased, only when the 
expansion in eq. (13) were to be truncated to some finite J < (L + l)^ to prevent the blowup of the eigenvalues A. 
Assuming for simplicity of the argument that J = the essence of the approach is now that the solution 



s(f) 



N 



(15) 



will be sparse (in achieving a bandwidth L using N Slepian instead of (L + 1)^ spherical harmonic expansion 
coefficients) yet good (in approximating the signal as best as possible in the mean squared sense within the 
region of interest R) and of geophysical utility (assuming we are dealing with spatially localized processes that 
are to be extracted, e.g., from global satellite measurements).^^ In light of the reasoning behind eqs (11)-(15), it 



is worth rereading the 1967 paper by W. M. Kaula,-*^^ which, written long before the advent of Slepian functions 
and the associated mathematical machinery, paved the way for many other studies.^' 

Bias and variance 

In concluding this section let us illustrate another welcome by-product of our methodology, by writing the mean 



squared error for the spherical harmonic solution (12) compared to the equivalent expression (14) for the solution 
in the Slepian basis. We do this as a function of the spatial coordinate, in the Slepian basis for both, and, for 
maximum clarity of the exposition, using the contrived case when both signal and noise should be bandlimited 
as well as white (both stipulations being technically impossible to satisfy simultaneously). In the former case. 



(L+l)^ (L+l)^ 

{e\r))=N ^ Xa[Xa^v{l-K)]-'gl{r)^v'S ^ (1 - A.)^[A. + r^(l - A.)]-^^^(f ), 



(16) 



while in the latter, we obtain 



N 



(L+l)2 



(17) 



a=l 



a>N 



1) basis functions are required to express the mean squared estimation error, whether in eq. (16) or in 



Ah (L 

eq. ( p!7| ). The ffist term in both expressions is the variance, which depends on the measurement noise. Without 
damping or truncation the variance grows without bounds. Damping and truncation alleviate this at the expense 
of added bias, which depends on the characteristics of the signal, as given by the second term. In contrast to 
T6l), however, the Slepian expression (17) has disentangled the contributions due to noise/ variance and 



eq 

signal/bias by projecting them onto the sparse set of well- localized and the remaining set of poorly localized 
Slepian functions, respectively. The estimation variance is felt via the basis functions a = 1 ^ N that are well 
concentrated inside the measurement area, and the effect of the bias is relegated to those a = N (L + l)^ 

functions that are confined to the region of missing data. 



When forming a solution to our problem in the Slepian basis by truncation according to eq. (15), changing 



the truncation level to values lower or higher than the Shannon number N amounts to navigating the trade-off 
space between variance, bias (or "resolution"), and sparsity in a manner that is captured with great clarity by 
eq. (17). We refer the reader elsewhere^ for more details. 



3. EXAMPLES AND APPLICATIONS 
Sparsity of errors in approximation 



Suppose that instead of eq. (11) we were to minimize 

/ {s — d)'^ dft = minimum. 



(18) 



which would lead to our forming the estimate from the partially observed data d directly as 



Sim = / dYimdVL. 

JR 



(19) 



While this might seem to be the natural approach in the absence of information outside the region of interest 



this would be equivalent to solving eq. (11) using a non-optimal damping constraint of = 1, as we have noted 
elsewhere,^ and eq. (12) furthermore shows this is so since Dim^i'm' Dim^i'm' = ^ii'^mm'- Nevertheless, we use 
this example because it is simple and informative, and it has been studied by others before. ^^'^^ If, for simplicity 
we assume that both signal and noise have a white power spectrum, the spectral error covariance matrix is 

{QmQ'm') = NDim^l'm' + SDim,l'm'^ (20) 

as can be derived by combining eq. (19) with eqs ([3| and ([7|-([8|. The errors are correlated, and are due to noise 
in the region R where data exist, and contaminated by signal from the missing region R = Q — R. It now makes 
immediate intuitive sense to desire a parameterization whose estimated coefficients have a diagonal covariance 
matrix. Let this parameterization be in terms of Slepian functions and thus this estimator 



Sa= dg^dn, (21) 
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which is to be contrasted with eq. ( p^O] ), and the associated error covariance will be diagonal and given by 

(6.6/3) = [NK + ^(1 - K)] Saf3^ (22) 

Eq. (22) follows from eq. (20) by the unitarity of the transform ([2| and the properties ([T]) and ([3|. 



An example of this behavior can be found in Figure [2j where the mean squared errors (ef^) and (e^) are 
plotted for the noiseless case, = 0, as a percentage of the signal strength, S. The geometry is that of a typical 
satellite survey characterized by a 'polar gap' R = Q — R consisting of a pair of axisymmetric caps at the North 
and South Pole, respectively, each one of various colatitudinal radii 6 = 0°, 3°, 5°, 7°, and 10°. The wedge shape 
of the errors in Figure [2]a-c? is bounded by the line \m\ = (/ + 1/2) sin©o. This relation substantiates earlier, 
heuristic, findings: it identifies Oq as the 'turning colatitude'^ separating the oscillatory from the evanescent 
parts of an asymptotic approximation of the Legendre functions at a given degree / and order m. Compared 



to the error structure of the spherical harmonic solution (19) shown in Figure pp-(i, the error structure of the 



Slepian basis solution (21) shown in Figure 2 is decidedly more sparse. 



Sparsity from geometry 

Geophysical signals that are regional in nature are sparse in the Slepian domain — provided the Slepian basis 
constructed is commensurate with the localization of the signal itself. To illustrate this we focus no longer on 
estimating the field but rather on the representation and approximation of the solution once it has been obtained. 
To this end we will drop the hats, omit the angular brackets on all symbols, and ignore observational noise. 

A bandlimited signal can be represented equally well in the spherical-harmonic as in any kind of Slepian 
basis of the same bandlimitation. When the signal is local and the chosen Slepian basis is similarly localized, 
approximating the function by the expansion truncated at the Shannon number will be advantageous 

L I (^+1)^ 

= ^ ^Irnyimi^) = ^ S^Qai^) - ^S^Qai^)^ (23) 

m= — l a=l a=l 

both because the quality of the approximation will be high in the region of interest and because the Shannon 
number of terms contributing to the reconstruction will be a significant savings over the full number of expansion 
coefficients. As per eq. (|6| this sparsity is thus mostly "geometric" in origin and the efficiency gains are dependent 
on the area of the region of interest expressed as a fraction of the area of the entire unit sphere. 

Using the double orthogonality definitions ([5| it is easy to derive from eq. ( 23 ) that the mean squared error 



(mse) over the region of interest of such an approximation will depend on the signal terms neglected in the 
expansion, and that the "i?-average mse" as a fraction of the i?-average mean signal strength should be given by 



/ e\r)dQ/ / s\r)dn= ^ ^A, / ^ ^A,. (24) 



The quahty of the regional approximation is high, because ~ 1 when a < N and Aq, ~ when a > N. 

Figure [3] illustrates this by expanding the regional "Bangui anomaly" in the lithospheric magnetic field in both 
the spherical-harmonic and Slepian bases. As can be seen in Figure |3]a- 6, the radial component of the Earth's 
main field, shown according to the POMME model^^ bandpass filtered to between degrees 17 and 72, contains 
some very prominent energy near Bangui, the capital of the Central African Republic. The origin of this anomaly 
remains debated. Windowing a spatial expansion of the field between degrees 17 and 36 with a Slepian window 
of L = 36 concentrated to a circular patch of radius © = 18° results in the rendition shown in Figure [3]c-g?. 
In the spherical harmonic domain, the windowed anomaly now contains energy between degrees / = — 72, as 
can be easily derived.^ Almost 80% of the spherical harmonic expansion coefficients are significant in that their 
absolute value rises above a threshold of one thousandth of their maximum absolute value. Switching bases and 
expanding the signal in the L = 72 Slepian basis localized to the same O = 18° region (i.e. those portrayed in 
Figure [T]) results in a very small number of significant expansion coefficients. Indeed, the partial reconstruction 
using the ffrst = 130 basis functions reveals that the anomaly is extremely well captured inside of the region 
of interest while the contribution to the signal is comprised of the energy of a mere 41 coefficients that are 
significant according to the same criterion, as shown in Figure [3]e-/. 

Sparsity from (geo) physics 

It has long been known that earthquakes perturb the terrestrial gravity field. Recently, thanks to the time- 
variable gravity measurements of the GRACE satellite pair,^^ the study of such changes has received renewed 
attention. ^^'^"^^^ Following seismological convention,^ let us denote the hypocentei of an earthquake, which is 
to be considered as a point source, as = {rg^Os^ (j)s)^ and let us vectorize the symmetric 'moment tensor' as 

M=[Mrr Mee M^^ Mro Me^]. (25) 

In a coordinate system = (r^O'^cj)') that is centered on the epicenter of the earthquake, the first-order Eule- 
rian gravitational potential perturbation in a spherically-symmetric non-rotating Earth is given by a sum over 
'spheroidal normal-mode' eigenfunctions ^P^(r) (which depend on the Earth model, in our case prem^'^), 

^^\r') = M- Y,T.n^r\Piir) i-^) nAi{r,,e',n (26) 

n=0 ^=0 ^ ^ 

where / is the usual spherical harmonic degree, the temporal frequency, and the vector excitation amplitude 

min(2,0 

nAl{rs,0\(j)') = ^ [nCirn{rs) COS mcj)' ^ nSlrn{rs)s^^rn(l)']Pirn{cOsO'), (27) 

m=0 

with m the usual spherical harmonic order, and with the associated Legendre function, P^^, evaluated at the epi- 
central angular distance 6'. The vector functions n^im and n^im are combinations of normal-mode displacement 
eigenfunctions and their radial derivatives,^' and need to be precomputed in the Earth model of choice. 

Solutions ( 26 ) for a variety of end- member earthquake sources are shown in Figure |4] As suggested by the 



summation limits of eq. (27), the patterns with which earthquakes perturb the Earth's gravity field have the 
symmetries of monopoles, dipoles, and quadrupoles.^' They should thus be eminently suitable to representation 
and amenable to analysis in the Slepian basis, in which, in other words, this type of geophysical signal is sparse. 

To the standard treatment using spherical Slepian functions we add one sophistication, namely rotation about 
their center of figure. As we have noted the Slepian basis set on circularly symmetric domains is degenerate and 
separable into 'colatitudinal' eigenfunctions which are solutions to a Sturm-Liouville equation that can be solved 
spectrally with great ease,^^ and order-dependent 'longitudinal' functions that control the axial symmetry. To 
compute the functions shown in Figure [l] we first determined the spherical harmonic expansion coefficients of the 
solutions to eq. (|3| centered on the North Pole, and subsequently rotated those to the colatitude and longitude 
of the desired cap center. Now we shall rotate the resulting functions over a third Euler angle about their center. 



We begin by a bait-and-switch for convenience — compared to Section [T] we will now work in the basis of 
complex surface spherical harmonics/'^ according to which the real Slepian and complex spherical harmonic 
expansion coefficients of the real signal s(f), and s/^, respectively, are related via the unitary transform 

L (L+l)' 

Sa = '^gllm^lm ^ud Si^ = ^ goclmScx, (28) 
Im cx=l 

where the Qaim form an (L + 1)^ x (L + 1)^ complex matrix. We shall here be concerned with obtaining the 
'localized' coefficients only when the spherical harmonic expansion coefficients sim are "known" , e.g. by being 
given in the form of so-called 'Level-2' GRACE data products. ^^'^^ Alternatives to this approach that can work 
directly from the raw data collected are discussed elsewhere.^' ^^'^^ 

Next, we adorn the Slepian basis functions, g^y^ their spherical harmonic expansion coefficients, gaimi and the 
Slepian expansion coefficients of the signal, s^^ by three upper indices, cj, 9 and describing, respectively, the 
rotation about their center of symmetry, < a; < 27r, and its colatitude, < ^ < tt, and longitude, < (/> < 27r. 
Thus, the "mother" Slepian basis set, concentrated over the circular region R centered on the North Pole, 
bandlimited to L, and of Shannon number N ^ consists of the functions ^^^^(f); they spawn a family of functions 
^^^'^(f), a = 1, . . . , (L + 1)^, centered at the geographical locations (^, (j)) on the unit sphere and rotated by u. 
For convenience we write no superscripts when {(j)^O^uo) = (0,0,0). The triad {(p^O^uo) contains the three Euler 
angles that generate a Slepian basis of arbitrary orientation anywhere on the sphere: first, by rotation over uo 
around the z axis, then by 6 about the original and finally by (j) around the original z axis again. The spherical 
harmonic coefficients of the rotated g"^^ and the mother set g^ are then related by^'^ 

9tTm^ E V'rrln.'i.4>,eM9o.l^-, (29) 

m' — — l 

where ^^^^^/(^, O^uj) is a (2/ + 1) x (2/ + 1) unitary transformation matrix of 'generalized' Legendre functions, ^^'^^ 

T>''2m'{<^^0,^) = e"""d«^,(^)e^™>, (30) 

duly noting that d^^^,(6) = d^^,^[—6). It is now convenient^^'^^ to decompose a general rotation to ((/), 6^u) into 
one over (0 — 7r/2, — 7r/2, 0) followed by another over (0, 7r/2, a; + 7r/2). In that case eq. (29) can be rewritten as 



The expansion of the signal s(y) into the rotated Slepian basis g'^ (r) is given, as in eqs (13) and (28), by 

.(f) = ^ st'-gf-{r) and = ^ ^^'^^/n^. (32) 



Combining eqs (31)-(32) and rearranging the summations we can write the 'fast Slepian expansion' as a discrete 
Fourier transforni7^'^~to be calculated on a grid of frequencies by EFT of the operator Z^^/^//, 



L 



JyQuj _ rpa -imuj -im' e -im" 4> ('X'X^^ 

— 2-^ 2-^ 2-^ ^mm'm"^ ^ ^ •> l^^Oaj 

m= — L m' = — L m" — — L 
L 

-'-mm'm" — ^mm' ^2 ) \2 J yalm"^lm- yC>C>D) 



^=max(|m| ,|m' | jlm" |) 



Figure [5] applies this procedure to the time- variable gravity field as seen by GRACE for the Slepian basis 
of colatitudinal radius 6 = 10° and bandwidth L = 60 centered at = 85° and = 95° and for varying 



rotations co. The gravitational-potential perturbation due to the 2004 Sumatra- Andaman earthquake is a clearly 
visible step in the time series of the expansion coefficients belonging to the a = 1 best concentrated m = ±1 
Slepian functions, expressed as height anomalies in m above the Earth's reference geoid. The orientation of the 
best-fitting fault plane can be derived from their relative contributions. 

Figure [6j finally, shows our best overall estimate of the geoid perturbation due to the Sumatra- Andaman 
earthquake. Figure [6]a-6 show renderings of the signal projected onto the best-concentrated m = and m = ±1 
components, respectively, of a Slepian basis with 6 = 10° and L = 60 and varying center locations. Each 
pixel in Figure |6]a and each arrow in Figure ^ plots the results from a different Slepian basis centered at the 
applicable locations (^, (p) shown. Figure [6]c shows the results from an inversion for the statistically significant 
"coseismic" step-changes of the geoid due to the earthquake in each of the well-concentrated basis functions 
of the single Slepian basis with O = 10° and L = 60 centered at = 85° and 0o = 95°, allowing also for 
exponential "postseismic" relaxation. ^^'^^ Figure [6]g? shows the signal predicted independently by expanding the 
geoid change obtained from a combined seismological-geodynamical forward model of the earthquake rupture 
and the associated coseismic perturbations in a simplified Cartesian Earth model. ^^'^^ More ample discussion 
of these models and their differences will take place in the geophysical literature. 
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Figure 1. Bandlimited eigenfunctions g{0,(j)) that are optimally concentrated within a circularly symmetric domain of 
colatitudinal radius G = 18° centered on Oq = 85° and 00 = 18°. The bandwidth is L = 72 and the rounded Shannon 
number N = 130. The circle denotes the cap boundary. Blue is positive and red is negative and the color axis is symmetric, 
but the sign is arbitrary; regions in which the absolute value is less than one hundredth of the maximum value on the 
sphere are left white. 
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Figure 2. Error covariance of the solutions to the geodetic estimation problem for scalar data observed in a region R in 
the presence of a 'polar gap', an axisymmetric pair of antipodal polar caps R = Q — R of radii G = 0°, 3°, 5°, 7°, and 10°, 
as shown, (a-d) Diagonal elements of the error covariance matrix, eq. (20), of the spherical harmonic solution (19). (e-h) 
Diagonal elements of the error covariance matrix, eq. (22), of the L = 90 Slepian basis solution (21). The ordinate is the 
sum of the rank a of the Slepian function within a sequence of single absolute order and this order, \m\. The lines at 
Z = I ml in panels a-h mark the area of the space of spherical harmonics at degrees / = 0, . . . , L, and orders |m|=0,...,/, 
while the lines at \m\ = (/ + 1/2) sin Go in panels a-d delineate the approximate influence zone of the errors. 
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Figure 3. Global and local representations of the lithospheric magnetic field in the spherical-harmonic and Slepian bases, 
(a) Map of the radial component of the internal magnetic field at the Earth's surface according to version 4.2s of the 
POMME model, bandpassed between spherical harmonic degrees 17 and 72, and (6) the spherical harmonic coefficients 
themselves. In total 5040 coefficients are needed to represent the global field, of which 5027 exceed a "1/1000" significance 
threshold of one thousandth of the maximum absolute value of all coefficients. Values below this relative threshold are 
left white in this and all other panels, (c) Map of what is known as the "Bangui anomaly", a highly localized feature in 
central Africa. The anomaly was obtained by multiplying the global field, low-passed to degree 36, by the Slepian function 
of bandwidth 36 that is best concentrated to the circular area of radius G = 18°, shown in blue. The resulting localized 
field is now bandlimited to degree 72, as can be seen in (d). Of the 5329 coefficients necessary to represent this anomaly, 
4181 exceed the "1/1000" threshold, (e) An approximation of the same anomaly using the N — 130 best-localized of the 
5329 Slepian functions concentrated to the region and bandlimited to degree 72. Of the 130 coefficients only 41 exceed 
the "1/1000" threshold, as shown in (/), which has a truncated ordinate. The approximation in the region of interest is 
beyond reproach and the representation by the Slepian expansion, compared to the spherical harmonics, is truly sparse. 



M = [1 -1 0]/V2 



M = [1 -1 0]/\/2 



20° 



10° 



-10° 





M = [0 1 0]/\/2 



M = [0 1 0]/\/2 



20° 



10° 



-10° 





1 0]/\/2 



20° 



10° 



-10° 




[0 -l]/\/2 




100° 110° 



80° 90° 100° 110° 



Figure 4. The patterns of gravitational potential perturbation owing to fictitious and idealized deviatoric double-couple 
point-source earthquakes occurring at 30 km depth underneath the island of Sumatra in the spherically symmetric Earth 
model prem,^^ calculated from normal-mode theory^ complete to degree L = 60 . Blue is positive and red is negative; the 
color axes are symmetric. The components of the moment tensor, eq. (25), of these end-member cases are indicated by the 



title. (a)-(6) 45°-dip thrust faults. {c)-{d) Vertical dip-slip faults. (e)-(/) Vertical strike-slip faults. All focal mechanisms 
were normalized to unit scalar moment magnitude. The relative proportions of the magnitudes of the response shown in 
each row are 100%, 54%, and 34%, respectively. 
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Figure 5. Orientation selectivity of the Slepian basis functions and the signal from the 12/26/2004 Sumatra- Andaman 
earthquake. The 'Level- 2' time- variable gravity spherical harmonic coefficients from the Gravity Recovery And Climate 



basis of colatitudinal radius G = 10° and of bandwidth L 
Oo = 85° and 00 = 95°. A linear trend, annual and semiannual variations were removed prior to display. 



Experiments*^ (GRACE) were transformed via eq. rt32j to the expansion coefficients in a circularly symmetric Slepian 

60, centered on the northwestern tip of the island of Sumatra, 

The monthly 

varying contributions from the m = ±1 best-concentrated (a = 1) basis functions are shown after rotation of the basis 
over the angles a; = 0°, —30° and —50°, respectively. The final rotation projects almost all of the energy of the signal onto 
a single component. By this measure, the best-fitting estimate of the overall strike of the earthquake is N40°W, which is 
in good agreement with independent seismological observations of this complex rupture. 
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Figure 6. The signal from the 12/26/2004 Sumatra- Andaman earthquake. The analysis was carried out using identical 
monthly solutions from GRACE as in Figure [s] and with the same Slepian basis functions, for which G = 10° and L — 60, 
but with their centers shifted as discussed below. The processing was identical in that a linear trend and (semi) annual 
variations were fitted and removed before inverting the resulting time series for the step-change due to the earthquake 
(the "coseismic signal"), (a) Magnitude of the zonal signal, i.e. the pointwise spatial expansion of the m — Slepian 
coefficients of the geoidal perturbation obtained after inversion of the time series for a step increase or decrease, using 
different Slepian basis sets whose centers (^, 0) coincide with the pixels being shown. (6) Directionality of the earthquake 
signal as measured by the pointwise expansion of the vectorial magnitude and direction of the m = zbl, a = 1 Slepian 
coefficients of the geoid change obtained by inversion of the GRACE time series. Every arrow drawn corresponds to the 
solution in a different Slepian basis centered at the current location, as is the case for each pixel in Figure [6]a. Blue and 
red indicate that the positive and negative parts, respectively, of the best-fitting rotated Slepian function are to be found 
to the north of the arrow, (c) The best estimate of the coseismic geoid change due to the Sumatra- Andaman earthquake 
obtained from inversion of the Slepian expansion of the GRACE time series for a step change at the known earthquake 
time, also allowing for an exponential relaxation of the perturbation (the "postseismic signal" ^^'^^). The spatial expansion 
is complete: the single Slepian basis used was centered on Oq = 85° and 00 = 95°. {d) Spatial rendition of the geoid 
change predicted independently, from seismic and geodynamical modeling of the earthquake. This forward model is 
based on several simplifying assumptions,^^ and more research is needed in order to compare Figures l6t-c? geophysically. 



